library(mgcv)
library(lme4)
library(reshape)
library(gamm4)
library(LMERConvenienceFunctions)
source("R-functions.R")

dat100mean<- read.table("spr-data.txt", head=TRUE)

dat100mean$ITEM_ID <- as.factor(ifelse(is.na(dat100mean$ITEM_NUMBER), NA, paste(substring(dat100mean$ITEM_TYPE,1,1),dat100mean$ITEM_NUMBER)))
summary(dat100mean)
areadat <- subset(dat100mean, !is.na(LEFT_AREA) & !is.na(RIGHT_AREA))

# calculate item order in experiment
avgtimes <- cast(melt(subset(areadat, !is.na(ITEM_ID))[,c(8,22,1)], id=1:2),
PARTICIPANT+ITEM_ID~variable, min, na.rm=TRUE) # item, time, participant
ca <- cast(melt(avgtimes, id=1:2), ITEM_ID~PARTICIPANT+variable, mean, na.rm=TRUE)
a =2
while(a<26){
ca[,a]<- rank(ca[,a])
a<- a+1
}
dab<- melt(ca)[,1:3];
names(dab)[2]<-"itemOrder"
d3 <- merge(areadat, dab, by=c("ITEM_ID","PARTICIPANT"), all=F, sort=F)
areadat <- subset(d3, !is.na(LEFT_AREA) & !is.na(RIGHT_AREA))

# throw out pupil dilations smaller than 2.5 SD avg for that participant
PARTICIPANT <- levels(areadat$PARTICIPANT)
vecdat100sds <- as.vector(tapply(areadat$LEFT_AREA,areadat$PARTICIPANT, sd))
vecareadats <- as.vector(tapply(areadat$LEFT_AREA,areadat$PARTICIPANT, mean))
ladf <-data.frame(vecareadats, vecdat100sds, PARTICIPANT)
d2 <-merge(areadat, ladf, by="PARTICIPANT", all=T, sort=F)
areadat$LEFT_AREA_n <- (d2$LEFT_AREA - d2$vecareadats) / (2.5*d2$vecdat100sds)
areadat$LEFT_AREA_n <- ifelse(areadat$LEFT_AREA_n < -1, NA, areadat$LEFT_AREA_n)

vecdat100sds <- as.vector(tapply(areadat$RIGHT_AREA,areadat$PARTICIPANT, sd))
vecareadats <- as.vector(tapply(areadat$RIGHT_AREA,areadat$PARTICIPANT, mean))
ladf <-data.frame(vecareadats, vecdat100sds, PARTICIPANT)
d2 <-merge(areadat, ladf, by="PARTICIPANT", all=T, sort = F)
areadat$RIGHT_AREA_n <- (d2$RIGHT_AREA - d2$vecareadats) / (2.5*d2$vecdat100sds)
areadat$RIGHT_AREA_n <- ifelse(areadat$RIGHT_AREA_n < -1, NA, areadat$RIGHT_AREA_n)

noblink <- subset(areadat, !is.na(RIGHT_AREA_n) & !is.na(LEFT_AREA_n))
data <- noblink

#calculate time per sentence, 0 at crit region onset
time <- subset(data, !is.na(INT_AREA))
time$CRIT_REG <- time$EXPINDEX - time$INT_AREA
intareatime <- subset(time, CRIT_REG==0)
reftime <- aggregate(intareatime$TIME, by=list(intareatime$PARTICIPANT,intareatime$ITEM_ID), FUN=min)
names(reftime) <-c("PARTICIPANT","ITEM_ID","REF_TIME")
data2 <- merge(data, reftime, by=c("PARTICIPANT","ITEM_ID"),sort=F)
rm(time) 
rm(intareatime)
rm(reftime)
data2$CRIT_REG <-data2$TIME-data2$REF_TIME

#create subsets
data2$ANSWER <- as.factor(data2$EXPKEYPRESS)

dataGEND <- subset(data2, ITEM_TYPE=="GB" | ITEM_TYPE=="GG")
dataGEND$RELWDINDEX <- dataGEND$EXPINDEX - dataGEND$INT_AREA

dataSEM <- subset(data2, ITEM_TYPE=="SB" | ITEM_TYPE=="SG")
dataSEM$SEMWDINDEX <- dataSEM$EXPINDEX - dataSEM$INT_AREA

dataRC <- subset(data2, ITEM_TYPE=="ra" | ITEM_TYPE=="rn")
dataRC$RELWDINDEX <- dataRC$EXPINDEX - dataRC$INT_AREA

###############

#### SPR GEND #####

pdf("SPR-GEND-icaraw.pdf", width=30, height=15, pointsize=40)
par(mfcol=c(1,1), lwd=4)
#### left eye ica events #####
src <- subset(dataGEND, ITEM_TYPE=="GG" & !is.na(LEFT_ICA_EVENT))
orc <- subset(dataGEND, ITEM_TYPE=="GB"& !is.na(LEFT_ICA_EVENT))
# src and orc are just variable names for the two conditions
stimevals <- tapply(src$CRIT_REG,src$CRIT_REG, mean)
srcvals <- tapply(atanh(src$LEFT_ICA_EVENT)*3,src$CRIT_REG, mean)
otimevals <- tapply(orc$CRIT_REG,orc$CRIT_REG, mean)
orcvals <- tapply(atanh(orc$LEFT_ICA_EVENT)*3,orc$CRIT_REG, mean)
plot(stimevals,srcvals, xlim=c(-2000,4000),type="l",ylim=c(2,4),# xlim=c(-500,1500),type="b",
#main="SPR GEND: Number of rapid dilations", 
ylab="number of rapid dilations within 500ms",
xlab="time in ms, 0 is onset of grammatical gender mismatch region")
rect(700,1,1250,26, col="lightgray", border="lightgray")
points(stimevals,srcvals,type="l")
points(otimevals,orcvals,col="red",type="l")
#legend("topright", c("correct", "incorrect"), col=c("black","red"), lty=1)
####  right eye ica events #####
src <- subset(dataGEND, ITEM_TYPE=="GG" & !is.na(RIGHT_ICA_EVENT))
orc <- subset(dataGEND, ITEM_TYPE=="GB"& !is.na(RIGHT_ICA_EVENT))
stimevals <- tapply(src$CRIT_REG,src$CRIT_REG, mean)
srcvals <- tapply(atanh(src$RIGHT_ICA_EVENT)*3,src$CRIT_REG, mean)
otimevals <- tapply(orc$CRIT_REG,orc$CRIT_REG, mean)
orcvals <- tapply(atanh(orc$RIGHT_ICA_EVENT)*3,orc$CRIT_REG, mean)
points(stimevals,srcvals,type="l", lty=3)
points(otimevals,orcvals, col="red",type="l", lty=3)
legend("topright", c("correct, left eye", "incorrect, left eye",
"correct, right eye", "incorrect, right eye"),
col=c("black","red","black","red"), lty=c(1,1,3,3),lwd=4)
dev.off()


a<- 700
b <- a+600
source("R-functions.R")
botheyes <- dataGEND[,c(1,2,4,5,7,8,12,21,22,23,27)]
css <- getAggregation(botheyes,a,b,"GG")
#scss <-subset(css,  ydiff<24)

#right eye
ml <- glmer(rawica ~ it + (1 +it| ITEM_ID) + (1+it| PARTICIPANT), data=subset(css, eye=="RIGHT_ICA_EVENT"), family=poisson)
summary(ml)
confint(ml, method="Wald")

#left eye
ml2 <- glmer(rawica ~ it + (1 +it| ITEM_ID) + (1+it| PARTICIPANT), data=subset(css, eye=="LEFT_ICA_EVENT"), family=poisson)
# left eye doesn't converge, therefore remove item slope for condition.
ml2 <- glmer(rawica ~ it + (1 | ITEM_ID) + (1+it| PARTICIPANT), data=subset(css, eye=="LEFT_ICA_EVENT"), family=poisson)
summary(ml2)
confint(ml2, method="Wald")


#both eyes:
# main interest:
m2b <- glmer(rawica ~ it+ (1+it| ITEM_ID) +(1+it| PARTICIPANT/eye), data=css, family=poisson )
summary(m2b)
confint(m2b, method="Wald")
# include item order:
m2bi <- glmer(rawica ~ it*scale(itemOrder)+ (1+it| ITEM_ID) +(1+it| PARTICIPANT/eye), data=css, family=poisson )
summary(m2bi)
confint(m2bi, method="Wald")
m2bia <- glmer(rawica ~ it+scale(itemOrder)+ (1+it| ITEM_ID) +(1+it| PARTICIPANT/eye), data=css, family=poisson )
summary(m2bia)
confint(m2bia, method="Wald")
anova(m2b,m2bia)
anova(m2bia,m2bi)
#baseline:
m2bbase <- glmer(rawica ~ 1+ (1+it| ITEM_ID) +(1+it| PARTICIPANT/eye), data=css, family=poisson )
anova(m2b,m2bbase)

# final selected model:
summary(m2b)
confint(m2b, method="Wald")

# fixation and saccade
m2bp <- gamm4(rawica ~ s(xpos) +it , random=~(1+it| ITEM_ID) +(1+it| PARTICIPANT/eye), data=css, family=poisson )
summary(m2bp$mer)


##### SPR SEM #####

pdf("SPR-SEM-icaraw.pdf", width=30, height=15, pointsize=40)
par(mfcol=c(1,1), lwd=4)
#### left eye ica events #####
src <- subset(dataSEM, ITEM_TYPE=="SG" & !is.na(LEFT_ICA_EVENT))
orc <- subset(dataSEM, ITEM_TYPE=="SB"& !is.na(LEFT_ICA_EVENT))
stimevals <- tapply(src$CRIT_REG,src$CRIT_REG, mean)
srcvals <- tapply(atanh(src$LEFT_ICA_EVENT)*3,src$CRIT_REG, mean)
otimevals <- tapply(orc$CRIT_REG,orc$CRIT_REG, mean)
orcvals <- tapply(atanh(orc$LEFT_ICA_EVENT)*3,orc$CRIT_REG, mean)
plot(stimevals,srcvals, xlim=c(-2000,4000),type="l",ylim=c(2,4),
#main="SPR SEM: Number of rapid dilations", 
ylab="number of rapid dilations per 100ms",
xlab="time in ms, 0 is onset of semantic mismatch region")
rect(750,1,1250,26, col="lightgray", border="lightgray")
points(stimevals,srcvals,type="l")
points(otimevals,orcvals,col="red",type="l")
legend("topright", c("correct, left eye", "incorrect, left eye",
"correct, right eye", "incorrect, right eye"),
col=c("black","red","black","red"), lty=c(1,1,3,3),lwd=2)
####  right eye ica events #####
src <- subset(dataSEM, ITEM_TYPE=="SG" & !is.na(RIGHT_ICA_EVENT))
orc <- subset(dataSEM, ITEM_TYPE=="SB"& !is.na(RIGHT_ICA_EVENT))
stimevals <- tapply(src$CRIT_REG,src$CRIT_REG, mean)
srcvals <- tapply(atanh(src$RIGHT_ICA_EVENT)*3,src$CRIT_REG, mean)
otimevals <- tapply(orc$CRIT_REG,orc$CRIT_REG, mean)
orcvals <- tapply(atanh(orc$RIGHT_ICA_EVENT)*3,orc$CRIT_REG, mean)
points(stimevals,srcvals,type="l", lty=3)
points(otimevals,orcvals, col="red",type="l", lty=3)
dev.off()

a<- 700
b <- a+600
source("R-functions.R")
botheyes <- dataSEM[,c(1,2,4,5,7,8,12,21,22,23,27)]
css <- getAggregation(botheyes,a,b,"SG")
#scss <-subset(css,  ydiff<24)

#left eye
mlleft <- glmer(rawica ~ it + (1 | ITEM_ID) + (1+it| PARTICIPANT), data=subset(css, eye=="LEFT_ICA_EVENT"), family=poisson)
summary(mlleft)
mllefta <- glmer(rawica ~ it*scale(itemOrder) + (1 | ITEM_ID) + (1+it| PARTICIPANT), data=subset(css, eye=="LEFT_ICA_EVENT"), family=poisson)
summary(mllefta)
anova(mlleft,mllefta)
confint(mlleft, method="Wald")
mlleftbase <-update(mlleft, .~.-it)
anova(mlleft,mlleftbase)

#right eye
mlright <- glmer(rawica ~ it + (1 | ITEM_ID) + (1+it| PARTICIPANT), data=subset(css, eye=="RIGHT_ICA_EVENT"), family=poisson)
summary(mlright)
mlrighta <- glmer(rawica ~ it*scale(itemOrder) + (1 | ITEM_ID) + (1+it| PARTICIPANT), data=subset(css, eye=="RIGHT_ICA_EVENT"), family=poisson)
summary(mlrighta)
anova(mlrighta,mlright)
confint(mlright, method="Wald")
mlrightbase <- update(mlright, .~.-it)
anova(mlrightbase,mlright)

#both eyes:
m2<-glmer(rawica ~ it *scale(itemOrder)+ (1+it| ITEM_ID) +(1+it| PARTICIPANT/eye), data=cssa,family=poisson)
summary(m2)
m2b<- update(m2,.~.-it:scale(itemOrder))
anova(m2,m2b)
summary(m2b)
m2c<-update(m2b,.~.-scale(itemOrder))
anova(m2b,m2c)
anova(m2,m2c)
confint(m2, method="Wald")
m2cbase <- update(m2c, .~.-it)
anova(m2c,m2cbase)
#final model:
m2a<-glmer(rawica ~ it *scale(itemOrder)+ (1+it| ITEM_ID) +(1+it*scale(itemOrder)| PARTICIPANT)+(1| PARTICIPANT:eye), data=css,family=poisson)
#final model:
summary(m2a)
confint(m2a,method="Wald")
anova(m2,m2a)

#comparison to model with fixation position:
m2ap<-glmer(rawica ~ xpos+ydiff+it *scale(itemOrder)+ (1+it| ITEM_ID) +(1+it| PARTICIPANT)+(1| PARTICIPANT:eye), data=css,family=poisson)
summary(m2ap)
m2apg<-gamm4(rawica ~ s(ydiff)+it *scale(itemOrder), random=~ (1+it| ITEM_ID) +(1+it| PARTICIPANT)+(1| PARTICIPANT:eye), data=css,family=poisson)
summary(m2apg$mer)

###### SPR RC #####

pdf("SPR-RC-icaraw.pdf", width=30, height=15, pointsize=40)
par(mfcol=c(1,1), lwd=4)
#### left eye ica events #####
src <- subset(dataRC, ITEM_TYPE=="rn" & !is.na(LEFT_ICA_EVENT))
orc <- subset(dataRC, ITEM_TYPE=="ra"& !is.na(LEFT_ICA_EVENT))
stimevals <- tapply(src$CRIT_REG,src$CRIT_REG, mean)
srcvals <- tapply(atanh(src$LEFT_ICA_EVENT)*3,src$CRIT_REG, mean)
otimevals <- tapply(orc$CRIT_REG,orc$CRIT_REG, mean)
orcvals <- tapply(atanh(orc$LEFT_ICA_EVENT)*3,orc$CRIT_REG, mean)
plot(stimevals,srcvals, xlim=c(-2000,4000),ylim=c(2,4),type="l",
#main="SPR RC: Number of rapid dilations",
 ylab="number of rapid dilations per 100ms",
xlab="time in ms, 0 is onset of disambiguating region")
rect(750,1,1250,6, col="lightgray", border="lightgray")
points(stimevals,srcvals, type="l")
points(otimevals,orcvals,col="red",type="l")
####  right eye ica events #####
src <- subset(dataRC, ITEM_TYPE=="rn" & !is.na(RIGHT_ICA_EVENT))
orc <- subset(dataRC, ITEM_TYPE=="ra"& !is.na(RIGHT_ICA_EVENT))
stimevals <- tapply(src$CRIT_REG,src$CRIT_REG, mean)
srcvals <- tapply(atanh(src$RIGHT_ICA_EVENT)*3,src$CRIT_REG, mean)
otimevals <- tapply(orc$CRIT_REG,orc$CRIT_REG, mean)
orcvals <- tapply(atanh(orc$RIGHT_ICA_EVENT)*3,orc$CRIT_REG, mean)
points(stimevals,srcvals, type="l",lty=3)
points(otimevals,orcvals,col="red",type="l",lty=3)
legend("topright", c("SRC, left eye", "ORC, left eye",
"SRC, right eye", "ORC, right eye"),
col=c("black","red","black","red"), lty=c(1,1,3,3),lwd=4)
dev.off()

source("R-functions.R")
botheyes <- dataRC[,c(1,2,4,5,7,8,12,21,22,23,27)]
a<- 700
b <- a+600
css <- getAggregation(botheyes,a,b,"rn")

# left eye
m1lo <- glmer(rawica ~ it*scale(itemOrder) + scale(ydiff) +(1 +it| ITEM_ID) + (1+it| PARTICIPANT), data=subset(css, eye=="LEFT_ICA_EVENT"), family=poisson)
summary(m1lo)# <- reported model
m1l <- glmer(rawica ~ it*scale(itemOrder) + (1 +it| ITEM_ID) + (1+it*scale(itemOrder)| PARTICIPANT), data=subset(css, eye=="LEFT_ICA_EVENT"), family=poisson)
# m1l doesn't converge.
m1l2 <- glmer(rawica ~ it + (1 +it| ITEM_ID) + (1+it| PARTICIPANT), data=subset(css, eye=="LEFT_ICA_EVENT"), family=poisson)
summary(m1l2)
m1lbase <- update(m1l, .~.-it)
anova(m1l, m1lbase)

#right eye
m1ra <- glmer(rawica ~ it*scale(itemOrder) + scale(ydiff)+ (1 +it| ITEM_ID) + (1+it| PARTICIPANT), data=subset(css, eye=="RIGHT_ICA_EVENT" ), family=poisson)
summary(m1ra)
m1r <- glmer(rawica ~ it*scale(itemOrder) + scale(ydiff)+ (1 +it| ITEM_ID) + (1+it*scale(itemOrder)| PARTICIPANT), data=subset(css, eye=="RIGHT_ICA_EVENT" ), family=poisson)
summary(m1r)
m1rbase <- update(m1r, .~.-it)
anova(m1r,m1rbase)
m1r <- gamm4(rawica ~ it*scale(itemOrder)+s(xpos)+s(ypos)+s(xdiff)+s(ydiff) , random=~ (1 +it| ITEM_ID) + (1+it| PARTICIPANT), data=subset(css, eye=="RIGHT_ICA_EVENT"), family=poisson)
summary(m1r$mer)

m1all <- glmer(rawica ~ it*scale(itemOrder)+ scale(ydiff)+(1 +it| ITEM_ID) + (1+it| PARTICIPANT)+ (1| PARTICIPANT:eye), data=subset(css), family=poisson)
summary(m1all)

#both eyes (first half)
m1first <- glmer(rawica ~ it+ scale(ydiff)+(1 +it| ITEM_ID) + (1+it| PARTICIPANT)+ (1| PARTICIPANT:eye), data=subset(css, itemOrder<50), family=poisson)
summary(m1first)
confint(m1first, method="Wald")


##############################
# DUAL TASK DATA
#############################
rm(list=ls())
library(mgcv)
library(lme4)
#library(lmerTest)
library(reshape)
library(gamm4)
source("R-functions.R")


dat100meano <- read.table("dual-task-data.txt", head=TRUE)

dat100meano$cond <- as.factor(ifelse(is.na(dat100meano$CONDITION), as.factor("-"), as.factor(dat100meano$CONDITION)))
levels(dat100meano$cond) <- levels(dat100meano$CONDITION)


dat100meano$phaselevel <- as.factor(dat100meano$phase)
data <- subset(dat100meano, !is.na(subject_name))
dat100meano <- data
starttimes <- as.vector(tapply(dat100meano$absolute_time, list(dat100meano$subject_name, dat100meano$phaselevel), min))
vecphaselevel <- as.vector(tapply(as.numeric(dat100meano$phaselevel), list(dat100meano$subject_name,dat100meano$phaselevel), mean))
vecpeople <- as.vector(tapply(as.numeric(dat100meano$subject_name), list(dat100meano$subject_name,dat100meano$phaselevel), mean))
subject_name <- levels(dat100meano$subject_name)[vecpeople]
phaselevel <- levels(dat100meano$phaselevel)[vecphaselevel]

phasestarttimes <- data.frame(starttimes, phaselevel, subject_name)

data <-merge(dat100meano, phasestarttimes, by=c("subject_name","phaselevel"), all=T, sort=F)
data$phasetimeinphase <- data$absolute_time - data$starttimes
dat100meano <- data
rm(data)

###calculate lateral target acceleration and lateral steering acceleration

veldat <-dat100meano[,c(5,8,13)]#expt2
veldat2 <- veldat
veldat2$absolute_time <- veldat2$absolute_time+100
accdat <- merge(veldat, veldat2, by=c("absolute_time"), sort=F)
accdat$target_acceleration <- accdat$lateral_target_velocity.x - accdat$lateral_target_velocity.y
accdat$steering_acceleration <- accdat$lateral_steering_velocity.x - accdat$lateral_steering_velocity.y
dat100mean2 <- merge(dat100meano, accdat[,c(1,6,7)], by=c("absolute_time"), sort=FALSE)
dat100mean <- dat100mean2
dat100mean$abs_target_acc <- abs(dat100mean$target_acceleration)
dat100mean$abs_steering_acc <- abs(dat100mean$steering_acceleration)
rm(dat100mean2)
rm(accdat)
rm(veldat2)
rm(veldat)


######################################################
### shift target bar velocity by 1.3 sec wrt. ICA
### shift steering bar velocity by 0.4 sec wrt. ICA
source("R-functions.R")
shift100data <- dat100mean[,c(1,8,35)]#expt2
shift100data$absolute_time <- shift100data$absolute_time+1300
names(shift100data)[2:3] <- c("SHIFT_LTV", "SHIFT_ALTA")

shift100data400 <- dat100mean[,c(1,8,13,35)]#expt2
shift100data400$absolute_time <- shift100data400$absolute_time+400
names(shift100data400)[2:4] <- c("SHIFT_LTV400", "SHIFT_LSV400", "SHIFT_ALTA400")

shift100data1600 <- dat100mean[,c(1,10)]#expt2
shift100data1600$absolute_time <- shift100data1600$absolute_time+1600
names(shift100data1600)[2] <- c("SHIFT_steerdev1600")

shifts1 <- merge(shift100data, shift100data1600, by="absolute_time", sort=F)
shifts <- merge(shifts1, shift100data400, by="absolute_time", sort=F)
dat100meano <- merge(dat100mean, shifts, by="absolute_time",sort=F)
rm(shift100data400)
rm(shift100data)
rm(shifts)
areadat <- subset(dat100meano, !is.na(LEFT_AREA) & !is.na(RIGHT_AREA))

vecphaselevel <- as.vector(tapply(as.numeric(areadat$phaselevel), list(areadat$subject_name,areadat$phaselevel), mean))
phaselevel <- levels(areadat$phaselevel)[vecphaselevel]
vecpeople <- as.vector(tapply(as.numeric(areadat$subject_name), list(areadat$subject_name,areadat$phaselevel), mean))
subject_name <- levels(areadat$subject_name)[vecpeople]
vecdat100sds <- as.vector(tapply(areadat$LEFT_AREA,list(areadat$subject_name, areadat$phaselevel), sd))
vecareadats <- as.vector(tapply(areadat$LEFT_AREA,list(areadat$subject_name, areadat$phaselevel), mean))
ladf <-data.frame(vecareadats, vecdat100sds, phaselevel, subject_name)

d2 <-merge(areadat, ladf, by=c("subject_name","phaselevel"), all=T, sort=F)
areadat$LEFT_AREA_n <- (d2$LEFT_AREA - d2$vecareadats) / (2.5*d2$vecdat100sds)
areadat$LEFT_AREA_n <- ifelse(areadat$LEFT_AREA_n < -1, NA, areadat$LEFT_AREA_n)

par(mfrow=c(4,6))
for (i in levels(areadat$subject_name)){
hist(subset(areadat, subject_name==i)$LEFT_AREA_n, breaks=100, main=i,xlim=c(-1.,1))
}

vecphaselevel2 <- as.vector(tapply(as.numeric(areadat$phaselevel), list(areadat$subject_name,areadat$phaselevel), mean))
phaselevel <- levels(areadat$phaselevel)[vecphaselevel]
vecpeople2 <- as.vector(tapply(as.numeric(areadat$subject_name), list(areadat$subject_name,areadat$phaselevel), mean))
subject_name <- levels(areadat$subject_name)[vecpeople]
vecdat100sdsr <- as.vector(tapply(areadat$RIGHT_AREA,list(areadat$subject_name, areadat$phaselevel), sd))
vecareadatsr <- as.vector(tapply(areadat$RIGHT_AREA,list(areadat$subject_name, areadat$phaselevel), mean))
ladf2 <-data.frame(vecareadatsr, vecdat100sdsr, phaselevel, subject_name)

d3 <-merge(areadat, ladf2, by=c("subject_name","phaselevel"), all=T, sort=F)
areadat$RIGHT_AREA_n <- (d3$RIGHT_AREA - d3$vecareadatsr) / (2.5*d3$vecdat100sdsr)
areadat$RIGHT_AREA_n <- ifelse(areadat$RIGHT_AREA_n < -1, NA, areadat$RIGHT_AREA_n)

noblink <- subset(areadat, !is.na(RIGHT_AREA_n) & !is.na(LEFT_AREA_n))

data2 <- dat100meano
data2noblink <- noblink

write.table(data2,"expt2data.txt", sep="\t")
write.table(data2noblink, "expt2data.noblink.txt", sep="\t")

rm(list=ls()[ls()!="dat100mean" & ls()!="getAggregationDual"])
# determine item order
noblinkexp2 <- read.table("expt2data.noblink.txt", head=TRUE)
noblinkexp2$time <- noblinkexp2$phasetimeinphase +( noblinkexp2$phaselevel * max(noblinkexp2$phasetimeinphase))
noblinkexp2$ITEM_ID <- as.factor(ifelse(noblinkexp2$ITEM_ID=="-", NA, paste(substring(noblinkexp2$CONDITION,1,1),noblinkexp2$ITEM_ID)))
avgtimes <- cast(melt(subset(noblinkexp2, !is.na(ITEM_ID))[,c(2,19,1)], id=1:2),
subject_name+ITEM_ID~variable, min, na.rm=TRUE) # item, time, participant
ca <- cast(melt(avgtimes, id=1:2), ITEM_ID~subject_name+variable, mean, na.rm=TRUE)
a =2
while(a<25){
ca[,a]<- rank(ca[,a])
a<- a+1
}
dab<- melt(ca)[,1:3];
names(dab)[2]<-"itemOrder"
d3 <- merge(noblinkexp2, dab, by=c("ITEM_ID","subject_name"), all=F, sort=F)
noblinkexp2 <- d3
rm(list=ls()[ls()!="dat100mean" & ls()!="noblinkexp2" & ls()!="getAggregationDual"])

grammardata <- subset(noblinkexp2, CONDITION=="gb" | CONDITION=="gg")
grammardata$CRIT_REG <- round(grammardata$TIME_SINCE_CRIT_ONSET / 100) * 100

semdata <- subset(noblinkexp2, CONDITION=="sb" | CONDITION=="sg")
semdata$CRIT_REG <- round(semdata$TIME_SINCE_CRIT_ONSET / 100) * 100

rcdatae2 <- subset(noblinkexp2, CONDITION=="rn" | CONDITION=="ra")
rcdatae2$CRIT_REG <- round(rcdatae2$TIME_SINCE_CRIT_ONSET / 100) * 100

drivingOnlye2 <- subset(noblinkexp2, CONDITION=="-" & is.na(TIME_SINCE_CRIT_ONSET))

source("R-functions.R")
#### DUAL GENDER ####

pdf("DUAL-GEND-icaraw.pdf", width=30, height=15, pointsize=40)
par(mfcol=c(1,1), lwd=4)
#### left eye ica events #####
src <- subset(grammardata, CONDITION=="gg" & !is.na(LEFT_ICA_EVENT))
orc <- subset(grammardata, CONDITION=="gb"& !is.na(LEFT_ICA_EVENT))
stimevals <- tapply(src$CRIT_REG,src$CRIT_REG, mean)
srcvals <- tapply(atanh(src$LEFT_ICA_EVENT)*3,src$CRIT_REG, mean)
otimevals <- tapply(orc$CRIT_REG,orc$CRIT_REG, mean)
orcvals <- tapply(atanh(orc$LEFT_ICA_EVENT)*3,orc$CRIT_REG, mean)
plot(stimevals,srcvals, xlim=c(-2000,4000),ylim=c(2,4),type="l",
#main="DUAL GRAMMAR: Number of rapid dilations", 
ylab="number of rapid dilations",
xlab="time in ms, 0 is onset of word with incorrect gender")
rect(750,1,1250,6, col="lightgray", border="lightgray")
points(stimevals,srcvals, type="l")
points(otimevals,orcvals,col="red",type="l")
####  right eye ica events #####
src <- subset(grammardata, CONDITION=="gg" & !is.na(RIGHT_ICA_EVENT))
orc <- subset(grammardata, CONDITION=="gb"& !is.na(RIGHT_ICA_EVENT))
stimevals <- tapply(src$CRIT_REG,src$CRIT_REG, mean)
srcvals <- tapply(atanh(src$RIGHT_ICA_EVENT)*3,src$CRIT_REG, mean)
otimevals <- tapply(orc$CRIT_REG,orc$CRIT_REG, mean)
orcvals <- tapply(atanh(orc$RIGHT_ICA_EVENT)*3,orc$CRIT_REG, mean)
points(stimevals,srcvals, type="l", lty=3)
points(otimevals,orcvals,col="red",type="l",lty=3)
legend("topright", c("correct, left eye", "incorrect, left eye",
"correct, right eye", "incorrect, right eye"),
col=c("black","red","black","red"), lty=c(1,1,3,3),lwd=4)
dev.off()

a<- 700
b<- a+600
botheyes <- grammardata[,c(2,1,30:33,20,16,18,50,51,45)]
css<-getAggregationDual(botheyes, a,b,"gg")


m1<-glmer(rawica ~ cond +scale(lsv_400)+ (1+cond| ITEM_ID) +(1+cond| subject_name)+(1| subject_name:eye), data=css, family=poisson)
summary(m1)
m1base <- update(m1, .~.-cond)
anova(m1,m1base)
m2<-glmer(rawica ~ cond +scale(lsv_400) + scale(xdiff) + scale(ydiff)+(1+cond| ITEM_ID) +(1+cond| subject_name)+(1| subject_name:eye), data=css, family=poisson)
summary(m2)
anova(m1,m2)
confint(m2, method="Wald")

#left eye
m2<-glmer(rawica ~ cond +scale(lsv_400)  + (1| ITEM_ID) +(1+cond| subject_name), data=subset(css, eye=="LEFT_ICA_EVENT"), family=poisson)
summary(m2)
m2base <- update(m2, .~.-cond)
anova(m2,m2base)

#right eye
m2<-glmer(rawica ~ cond+scale(itemOrder) +scale(lsv_400)  + (1+cond| ITEM_ID) +(1+cond| subject_name), data=subset(css, eye=="RIGHT_ICA_EVENT"), family=poisson)
summary(m2)
m2base <- update(m2, .~.-cond)
anova(m2,m2base)


rm(list=c("a", "b","bl","botheyes","css","m1", "m1b","m2","m2b","mbotheyes","mf","orc","src","orcvals","otimevals","src","srcvals","stimevals"))

#### DUAL SEM #####

pdf("DUAL-SEM-icaraw.pdf", width=30, height=15, pointsize=40)
par(mfcol=c(1,1), lwd=4)
#### left eye ica events #####
src <- subset(semdata, CONDITION=="sg" & !is.na(LEFT_ICA_EVENT)& ITEM_ID != "s 26"& ITEM_ID != "s 6"& ITEM_ID != "s 15")
orc <- subset(semdata, CONDITION=="sb"& !is.na(LEFT_ICA_EVENT)& ITEM_ID != "s 26"& ITEM_ID != "s 6"& ITEM_ID != "s 15")
stimevals <- tapply(src$CRIT_REG,src$CRIT_REG, mean)
srcvals <- tapply(atanh(src$LEFT_ICA_EVENT)*3,src$CRIT_REG, mean)
otimevals <- tapply(orc$CRIT_REG,orc$CRIT_REG, mean)
orcvals <- tapply(atanh(orc$LEFT_ICA_EVENT)*3,orc$CRIT_REG, mean)
plot(stimevals,srcvals, xlim=c(-2000,4000),ylim=c(2,4),type="l",
#main="DUAL SEM: Number of rapid dilations", 
ylab="number of rapid dilations",
xlab="time in ms, 0 onset of the semantically anomalous word")
rect(750,1,1250,6, col="lightgray", border="lightgray")
rect(450,1,950,6, border="black", lty=2)
points(stimevals,srcvals, type="l")
points(otimevals,orcvals,col="red",type="l")
####  right eye ica events #####
src <- subset(semdata, CONDITION=="sg" & !is.na(RIGHT_ICA_EVENT)& ITEM_ID != "s 26"& ITEM_ID != "s 6"& ITEM_ID != "s 15")
orc <- subset(semdata, CONDITION=="sb"& !is.na(RIGHT_ICA_EVENT)& ITEM_ID != "s 26"& ITEM_ID != "s 6"& ITEM_ID != "s 15")
stimevals <- tapply(src$CRIT_REG,src$CRIT_REG, mean)
srcvals <- tapply(atanh(src$RIGHT_ICA_EVENT)*3,src$CRIT_REG, mean)
otimevals <- tapply(orc$CRIT_REG,orc$CRIT_REG, mean)
orcvals <- tapply(atanh(orc$RIGHT_ICA_EVENT)*3,orc$CRIT_REG, mean)
points(stimevals,srcvals, type="l",lty=3)
points(otimevals,orcvals,col="red",type="l",lty=3)
legend("topright", c("correct, left eye", "incorrect, left eye",
"correct, right eye", "incorrect, right eye"),
col=c("black","red","black","red"), lty=c(1,1,3,3),lwd=4)
dev.off()

a<- 400
b<- a+600
botheyes <- semdata[,c(2,1,30:33,20,16,18,50,51,45)]
css<-getAggregationDual(botheyes, a,b,"sg")
css <- subset(css, ITEM_ID != "s 26"& ITEM_ID != "s 6"& ITEM_ID != "s 15")

m1<-glmer(rawica ~ cond+scale(ydiff) + (1+cond| ITEM_ID) +(1+scale(ydiff) +cond| subject_name)+(1| subject_name:eye), data=subset(css), family=poisson)
summary(m1)
confint(m1, method="Wald")
m1base <- update(m1, .~.-cond)
anova(m1,m1base)


#left eye
m2<-glmer(rawica ~ cond +scale(lsv_400) + scale(xdiff) + scale(ydiff)+(1+cond| ITEM_ID) +(1+cond| subject_name), data=subset(css, eye="LEFT_ICA_EVENT"), family=poisson)
summary(m2)
confint(m2, method="Wald")

#left eye
m2<-glmer(rawica ~ cond +scale(lsv_400)  + (1+cond| ITEM_ID) +(1+cond| subject_name), data=subset(css, eye=="LEFT_ICA_EVENT"), family=poisson)
summary(m2)
m2base <- update(m2, .~.-cond)
anova(m2,m2base)

#right eye
m2<-glmer(rawica ~ cond +scale(lsv_400) + (1+cond| ITEM_ID) +(1+cond| subject_name), data=subset(css, eye=="RIGHT_ICA_EVENT"), family=poisson)
summary(m2)
m2base <- update(m2, .~.-cond)
anova(m2,m2base)


#residdata<- subset(mbotheyes, CRIT_REG > a & CRIT_REG < b & !is.na(rawica) &RIGHT_X > 0 & RIGHT_Y > 0 & RIGHT_X < 1960 & RIGHT_Y <1200 )



rm(list=c("a", "b","botheyes","css","m1","m2","orcvals","otimevals","src","srcvals","stimevals"))

#### DUAL RC ####

pdf("DUAL-RC-icaraw.pdf", width=30, height=15, pointsize=40)
par(mfcol=c(1,1), lwd=4)
#### left eye ica events #####
src <- subset(rcdatae2, CONDITION=="rn" & !is.na(LEFT_ICA_EVENT))
orc <- subset(rcdatae2, CONDITION=="ra"& !is.na(LEFT_ICA_EVENT))
stimevals <- tapply(src$CRIT_REG,src$CRIT_REG, mean)
srcvals <- tapply(atanh(src$LEFT_ICA_EVENT)*3,src$CRIT_REG, mean)
otimevals <- tapply(orc$CRIT_REG,orc$CRIT_REG, mean)
orcvals <- tapply(atanh(orc$LEFT_ICA_EVENT)*3,orc$CRIT_REG, mean)
plot(stimevals,srcvals, xlim=c(-2000,4000),ylim=c(2,4),type="l",
#main="DUAL RC: Number of rapid dilations", 
ylab="number of rapid dilations",
xlab="time in ms, 0 is onset of disambiguating region")
rect(750,1,1250,6, col="lightgray", border="lightgray")
points(stimevals,srcvals, type="l")
points(otimevals,orcvals,col="red",type="l")
####  right eye ica events #####
src <- subset(rcdatae2, CONDITION=="rn" & !is.na(RIGHT_ICA_EVENT))
orc <- subset(rcdatae2, CONDITION=="ra"& !is.na(RIGHT_ICA_EVENT))
stimevals <- tapply(src$CRIT_REG,src$CRIT_REG, mean)
srcvals <- tapply(atanh(src$RIGHT_ICA_EVENT)*3,src$CRIT_REG, mean)
otimevals <- tapply(orc$CRIT_REG,orc$CRIT_REG, mean)
orcvals <- tapply(atanh(orc$RIGHT_ICA_EVENT)*3,orc$CRIT_REG, mean)
points(stimevals,srcvals, type="l", lty=3)
points(otimevals,orcvals,col="red",type="l", lty=3)
legend("topright", c("SRC, left eye", "ORC, left eye","SRC, right eye", "ORC, right eye"), col=c("black","red","black","red"), lty=c(1,1,3,3),lwd=4)
dev.off()

a<- 700
b<- a+600
botheyes <- rcdatae2[,c(2,1,30:33,20,16,18,50,51,45)]
css<-getAggregationDual(botheyes, a,b,"rn")

m1<-glmer(rawica ~ cond +scale(lsv_400)+(1+cond| ITEM_ID) +(1 +cond+scale(lsv_400)| subject_name)+(1| subject_name:eye), data=subset(css), family=poisson)
summary(m1)
confint(m1, method="Wald")
m1base <- update(m1, .~.-cond)
anova(m1,m1base)


#left eye
m2<-glmer(rawica ~ cond +scale(lsv_400)+(1+cond| ITEM_ID) +(1 +cond| subject_name), data=subset(css, eye="LEFT_ICA_EVENT"), family=poisson)
summary(m2)
confint(m2, method="Wald")
m2base <- update(m2, .~.-cond)
anova(m2,m2base)

#right eye
m2<-glmer(rawica ~ cond +scale(lsv_400)+(1+cond| ITEM_ID) +(1 +cond| subject_name), data=subset(css, eye=="RIGHT_ICA_EVENT"), family=poisson)
summary(m2)
m2base <- update(m2, .~.-cond)
anova(m2,m2base)


############
rcdatae2$leftrawica <- atanh(rcdatae2$LEFT_ICA_EVENT)* 3
rcdatae2$rightrawica <- atanh(rcdatae2$RIGHT_ICA_EVENT)* 3
botheyes <- rcdatae2[,c(2,1,20,45,50,51,52,53,30:33)]
mbotheyes <- melt(botheyes, id= c(1:6,9:12) )
names(mbotheyes)[12]<-"rawica"
names(mbotheyes)[11]<-"eye"
a<- 700
b<- a+600
residdata<-subset(mbotheyes, CRIT_REG > a & CRIT_REG < b & !is.na(rawica)&RIGHT_X > 0 & RIGHT_Y > 0 & RIGHT_X < 1960 & RIGHT_Y <1200)
bl<-lm(rawica~SHIFT_LSV400 ,data= residdata)
residica <- resid(bl)
summary(bl)
confint(bl, method="Wald")
residdata$residica <- residica
residcomb <- residdata[,c(1:3,5,6,11,13)]
residcomb2 <- melt(residcomb, id= 1:6)  #c("PARTICIPANT", "ITEM_TYPE","ITEM_ID", "CRIT_REG"))
css<-cast(residcomb2, subject_name + CONDITION + ITEM_ID +itemOrder+eye~ variable,
mean)
css$cond <- ifelse(css$CONDITION=="rn", 0, 1)
m1<-lmer(residica ~ cond *itemOrder+ (1| ITEM_ID) +(1| subject_name), data=css)
m2<-fitLMER.fnc(model=m1, ran.effects=c("(1| eye)","(0+cond|ITEM_ID)","(0+cond| subject_name)","(0+cond| eye)"),backfit.on="F",log.file=FALSE)
summary(m2)


mf1 <- lmer(residica ~ cond + (1+cond | ITEM_ID) +
(1+cond| subject_name)+(1+cond| eye), data=css)
summary(mf1)
confint(mf1,method="Wald")


##############################
## effect of driving task ######
##############################
rm(list=ls()[ls()!="dat100mean" ])

par(mfrow=c(4,6))
a<-1
while(a<24){
p <- levels(dat100mean$subject_name)[a]
sldat<-left_ICA_diff <- atanh(subset(dat100mean, subject_name==p & 
!is.na(lateral_target_velocity) & CONDITION=="-" & difficulty=="d"
& !is.na(LEFT_ICA_EVENT) & !is.na(RIGHT_ICA_EVENT))$LEFT_ICA_EVENT)*3
sldat2<-subset(dat100mean, subject_name==p & 
!is.na(lateral_target_velocity) & CONDITION=="-" & difficulty=="d"
& !is.na(LEFT_ICA_EVENT) & !is.na(RIGHT_ICA_EVENT))$lateral_target_velocity
ccf(sldat,sldat2, lag.max=50, main=a,xlim=c(-30,30))
a<-a+1;
}

dev.new()
ddd<-subset(dat100mean, !is.na(lateral_target_velocity) &
CONDITION=="-" & difficulty=="d" & !is.na(LEFT_ICA_EVENT) &
!is.na(RIGHT_ICA_EVENT))
subj_perf<- tapply(ddd$steering_deviation, ddd$subject_name, mean)
names(subj_perf)<- 1:23
barplot(subj_perf)

##################
## VISUAL WORLD ##
##################
rm(list=ls())
data <- read.table("visualWorld-data.txt", head=TRUE)
source("R-functions.R")
data$item <- as.factor(data$item)
data <- subset(data, condition != "14aiia" )# excluded before analysis
#due to error in item
data$partnum <- as.numeric(substr(data$participant,0,2))
areadat <- subset(data, !is.na(LEFT_AREA) & !is.na(RIGHT_AREA))

PARTICIPANT <- levels(areadat$PARTICIPANT)
vecdat100sds <- as.vector(tapply(areadat$LEFT_AREA,areadat$PARTICIPANT, sd))
vecareadats <- as.vector(tapply(areadat$LEFT_AREA,areadat$PARTICIPANT, mean))
ladf <-data.frame(vecareadats, vecdat100sds, PARTICIPANT)
d2 <-merge(areadat, ladf, by="PARTICIPANT", all=T, sort=F)
areadat$LEFT_AREA_n <- (d2$LEFT_AREA - d2$vecareadats) / (2.5*d2$vecdat100sds)
areadat$LEFT_AREA_n <- ifelse(areadat$LEFT_AREA_n < -1, NA, areadat$LEFT_AREA_n)

vecdat100sds <- as.vector(tapply(areadat$RIGHT_AREA,areadat$PARTICIPANT, sd))
vecareadats <- as.vector(tapply(areadat$RIGHT_AREA,areadat$PARTICIPANT, mean))
ladf <-data.frame(vecareadats, vecdat100sds, PARTICIPANT)
d2 <-merge(areadat, ladf, by="PARTICIPANT", all=T, sort = F)
areadat$RIGHT_AREA_n <- (d2$RIGHT_AREA - d2$vecareadats) / (2.5*d2$vecdat100sds)
areadat$RIGHT_AREA_n <- ifelse(areadat$RIGHT_AREA_n < -1, NA, areadat$RIGHT_AREA_n)

noblink <- subset(areadat, !is.na(RIGHT_AREA_n) & !is.na(LEFT_AREA_n))
data <- noblink
avgtimes <- cast(melt(subset(data, !is.na(item))[,c(1,6,7)], id=2:3),
participant+item~variable, min, na.rm=TRUE)
ca <- cast(melt(avgtimes, id=1:2), item~participant+variable, mean, na.rm=TRUE)
a =2
while(a<25){
ca[,a]<- rank(ca[,a])
a<- a+1
}
dab<- melt(ca)[,1:3];
names(dab)[2]<-"itemOrder"
d3 <- merge(data, dab, by=c("item","participant"), all=T, sort=F)
data<-d3
conntimes <- read.table("cond-connectorstart.txt")
names(conntimes)<- c("condition", "CRIT_ONSET")
d3<-merge(data, conntimes, by=c("condition"), all=T, sort=F)
data<-d3
data$CRIT_REG <- data$timediff - data$CRIT_ONSET
data$rCRIT_REG <- round(data$CRIT_REG, digits=1)
early <- subset(data, LEFT_X > 0 & LEFT_X < 1920 & LEFT_Y > 0 & LEFT_Y < 1200)
data <- early[,c(3,2,10,11,13,14,8,18,19,25,28)]

pdf("VW-icaraw.pdf", width=30, height=15, pointsize=40)
par(mfcol=c(1,1), lwd=4)
#### left eye ica events #####
src <- subset(data, concessive==FALSE & !is.na(LEFT_ICA_EVENT) & itemOrder<11)
orc <- subset(data, concessive==TRUE & !is.na(LEFT_ICA_EVENT) & itemOrder<11)
stimevals <- tapply(src$rCRIT_REG,src$rCRIT_REG, mean)
srcvals <- tapply(atanh(src$LEFT_ICA_EVENT)*3,src$rCRIT_REG, mean)
otimevals <- tapply(orc$rCRIT_REG,orc$rCRIT_REG, mean)
orcvals <- tapply(atanh(orc$LEFT_ICA_EVENT)*3,orc$rCRIT_REG, mean)
plot(stimevals,srcvals, ylim=c(2,4),type="l", xlim=c(-2,4),
main="Visual World Connective: Number of rapid dilations", ylab="number of rapid dilations",
xlab="time in ms, 0 is onset of connective")
rect(.750,1,1.250,6, col="lightgray", border="lightgray")
points(stimevals,srcvals, type="l")
points(otimevals,orcvals,col="red",type="l")
####  right eye ica events #####
src <- subset(data, concessive==FALSE & !is.na(RIGHT_ICA_EVENT)& itemOrder<11)
orc <- subset(data, concessive==TRUE& !is.na(RIGHT_ICA_EVENT)& itemOrder<11)
stimevals <- tapply(src$rCRIT_REG,src$rCRIT_REG, mean)
srcvals <- tapply(atanh(src$RIGHT_ICA_EVENT)*3,src$rCRIT_REG, mean)
otimevals <- tapply(orc$rCRIT_REG,orc$rCRIT_REG, mean)
orcvals <- tapply(atanh(orc$RIGHT_ICA_EVENT)*3,orc$rCRIT_REG, mean)
points(stimevals,srcvals, type="l",lty=3)
points(otimevals,orcvals,col="red",type="l",lty=3)
legend("topright", c("correct, left eye", "incorrect, left eye",
"correct, right eye", "incorrect, right eye"),cex=1,
col=c("black","red","black","red"), lty=c(1,1,3,3),lwd=4)
dev.off()

pdf("VW-3-icaraw.pdf", width=30, height=25, pointsize=40)
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
par(lwd=4)
#### left eye ica events #####
src <- subset(data, concessive==FALSE & !is.na(LEFT_ICA_EVENT))
orc <- subset(data, concessive==TRUE & !is.na(LEFT_ICA_EVENT))
stimevals <- tapply(src$rCRIT_REG,src$rCRIT_REG, mean)
srcvals <- tapply(atanh(src$LEFT_ICA_EVENT)*3,src$rCRIT_REG, mean)
otimevals <- tapply(orc$rCRIT_REG,orc$rCRIT_REG, mean)
orcvals <- tapply(atanh(orc$LEFT_ICA_EVENT)*3,orc$rCRIT_REG, mean)
plot(stimevals,srcvals, ylim=c(2,4),type="l", xlim=c(-2,4),
main="(a) Concessive vs. causal connective in VW study", ylab="number of rapid dilations",
xlab="time in ms, 0 is onset of connective")
rect(.750,1,1.250,6, col="lightgray", border="lightgray")
points(stimevals,srcvals, type="l")
points(otimevals,orcvals,col="red",type="l")
####  right eye ica events #####
src <- subset(data, concessive==FALSE & !is.na(RIGHT_ICA_EVENT))
orc <- subset(data, concessive==TRUE& !is.na(RIGHT_ICA_EVENT))
stimevals <- tapply(src$rCRIT_REG,src$rCRIT_REG, mean)
srcvals <- tapply(atanh(src$RIGHT_ICA_EVENT)*3,src$rCRIT_REG, mean)
otimevals <- tapply(orc$rCRIT_REG,orc$rCRIT_REG, mean)
orcvals <- tapply(atanh(orc$RIGHT_ICA_EVENT)*3,orc$rCRIT_REG, mean)
points(stimevals,srcvals, type="l",lty=3)
points(otimevals,orcvals,col="red",type="l",lty=3)
legend("topright", c("correct, left eye", "incorrect, left eye",
"correct, right eye", "incorrect, right eye"),cex=1,
col=c("black","red","black","red"), lty=c(1,1,3,3),lwd=4)
#### left eye ica events #####
src <- subset(data, concessive==FALSE & !is.na(LEFT_ICA_EVENT) & itemOrder<11)
orc <- subset(data, concessive==TRUE & !is.na(LEFT_ICA_EVENT) & itemOrder<11)
stimevals <- tapply(src$rCRIT_REG,src$rCRIT_REG, mean)
srcvals <- tapply(atanh(src$LEFT_ICA_EVENT)*3,src$rCRIT_REG, mean)
otimevals <- tapply(orc$rCRIT_REG,orc$rCRIT_REG, mean)
orcvals <- tapply(atanh(orc$LEFT_ICA_EVENT)*3,orc$rCRIT_REG, mean)
plot(stimevals,srcvals, ylim=c(2,4),type="l", xlim=c(-1,2),
main="(b) first half of expt", ylab="number of rapid dilations",
xlab="time in ms, 0 is onset of connective")
rect(.750,1,1.250,6, col="lightgray", border="lightgray")
points(stimevals,srcvals, type="l")
points(otimevals,orcvals,col="red",type="l")
####  right eye ica events #####
src <- subset(data, concessive==FALSE & !is.na(RIGHT_ICA_EVENT)& itemOrder<=10)
orc <- subset(data, concessive==TRUE& !is.na(RIGHT_ICA_EVENT)& itemOrder<=10)
stimevals <- tapply(src$rCRIT_REG,src$rCRIT_REG, mean)
srcvals <- tapply(atanh(src$RIGHT_ICA_EVENT)*3,src$rCRIT_REG, mean)
otimevals <- tapply(orc$rCRIT_REG,orc$rCRIT_REG, mean)
orcvals <- tapply(atanh(orc$RIGHT_ICA_EVENT)*3,orc$rCRIT_REG, mean)
points(stimevals,srcvals, type="l",lty=3)
points(otimevals,orcvals,col="red",type="l",lty=3)
#### left eye ica events #####
src <- subset(data, concessive==FALSE & !is.na(LEFT_ICA_EVENT) & itemOrder>10)
orc <- subset(data, concessive==TRUE & !is.na(LEFT_ICA_EVENT) & itemOrder>10)
stimevals <- tapply(src$rCRIT_REG,src$rCRIT_REG, mean)
srcvals <- tapply(atanh(src$LEFT_ICA_EVENT)*3,src$rCRIT_REG, mean)
otimevals <- tapply(orc$rCRIT_REG,orc$rCRIT_REG, mean)
orcvals <- tapply(atanh(orc$LEFT_ICA_EVENT)*3,orc$rCRIT_REG, mean)
plot(stimevals,srcvals, ylim=c(2,4),type="l", xlim=c(-1,2),
main="(c) second half of expt", ylab="number of rapid dilations",
xlab="time in ms, 0 is onset of connective")
rect(.750,1,1.250,6, col="lightgray", border="lightgray")
points(stimevals,srcvals, type="l")
points(otimevals,orcvals,col="red",type="l")
####  right eye ica events #####
src <- subset(data, concessive==FALSE & !is.na(RIGHT_ICA_EVENT)& itemOrder>10)
orc <- subset(data, concessive==TRUE& !is.na(RIGHT_ICA_EVENT)& itemOrder>10)
stimevals <- tapply(src$rCRIT_REG,src$rCRIT_REG, mean)
srcvals <- tapply(atanh(src$RIGHT_ICA_EVENT)*3,src$rCRIT_REG, mean)
otimevals <- tapply(orc$rCRIT_REG,orc$rCRIT_REG, mean)
orcvals <- tapply(atanh(orc$RIGHT_ICA_EVENT)*3,orc$rCRIT_REG, mean)
points(stimevals,srcvals, type="l",lty=3)
points(otimevals,orcvals,col="red",type="l",lty=3)
dev.off()

source("R-functions.R")
a<- .7
b<- a+.600
css <- getAggregationVW(data, a, b)

m1<-glmer(rawica ~ concessive*scale(itemOrder) +scale(ydiff)+scale(ypos)+(1 +concessive | item) +(1 + concessive*scale(itemOrder)| participant)+(1 | participant:eye), data=subset(css), family=poisson)
summary(m1) 
m2<-update(m1, .~.-concessive)
anova(m1,m2)

m1h<-glmer(rawica ~ concessive +scale(ydiff)+scale(ypos)+(1+ concessive | item) +(1+concessive | participant)+(1 | participant:eye), data=subset(css, itemOrder<11), family=poisson)
summary(m1h) 
confint(m1h, method="Wald")
m2h<-update(m1h, .~.-concessive)
anova(m1h,m2h)


m1r<-glmer(rawica ~ concessive+scale(ypos)+scale(ydiff)+(1+concessive  | item) +(1 + concessive| participant), data=subset(css, eye=="RIGHT_ICA_EVENT"& itemOrder<11), family=poisson)
summary(m1r) 
confint(m1r, method="Wald")

print(dotplot(ranef(m1,condVar=TRUE), scales = list(x = list(relation = 'free')))[["participant"]])
print(dotplot(ranef(m1,condVar=TRUE), scales = list(x = list(relation = 'free')))[["item"]])

m1 <-gamm4(rawica ~ concessive +s(xpos) + s(ypos) , random=~(1+ concessive | item) +(1+concessive|participant/eye), data=css, family = poisson)
summary(m1$mer)
plot(m1$gam) 
confint(m1, method="Wald")
m1base <- update(m1, .~.-cond)
anova(m1,m1base)


#left eye
m2<-glmer(rawica ~ cond +scale(lsv_400)+(1+cond| ITEM_ID) +(1 +cond| subject_name), data=subset(css, eye="LEFT_ICA_EVENT"), family=poisson)
summary(m2)
confint(m2, method="Wald")
m2base <- update(m2, .~.-cond)
anova(m2,m2base)

#right eye
m2<-glmer(rawica ~ cond +scale(lsv_400)+(1+cond| ITEM_ID) +(1 +cond| subject_name), data=subset(css, eye=="RIGHT_ICA_EVENT"), family=poisson)
summary(m2)
m2base <- update(m2, .~.-cond)
anova(m2,m2base)


